Spatially heterogeneous dynamics and dynamic facilitation in a model of viscous silica 
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Performing molecular dynamics simulations, we find that the structural relaxation dynamics of 
viscous silica, the prototype of a strong glass former, are spatially heterogeneous and cannot be 
understood as a statistical bond breaking process. Further, we show that high particle mobility 
predominantly propagates continuously through the melt, supporting the concept of dynamic facil- 
itation emphasized in recent theoretical work. 
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One of the most challenging problems of condensed 
matter physics is still the understanding of the tremen- 
dous slowing down of molecular dynamics in supercooled 
liquids approaching the glass transition temperature T g , 
which is not accompanied by a substantial change of the 
structure. In the past, several theories have been in- 
voked to rationalize this behavior. The mode coupling 
theory 1] describes many experimental findings at high 
temperatures, but it predicts a power-law divergence of 
relaxation times at a critical temperature Tmct signifi- 
cantly above T g . The Adam-Gibbs approach |2j sees the 
glass transition phenomenon as resulting from an increase 
of cooperativity of molecular dynamics when the temper- 
ature is decreased. However, the underlying microscopic 
picture is not well established. 

Very recently, Garrahan and Chandler (GC) 0, Q 
introduced a microscopic model of supercooled liquids, 
which is based on three central ideas: (i) Particle mo- 
bility is sparse and dynamics are spatially heterogeneous 
at times intermediate between ballistic and diffusive mo- 
tion, (ii) Particle mobility is the result of dynamic fa- 
cilitation, i.e., mobile particles assist their neighbors to 
become mobile, (hi) Mobility propagation carries a direc- 
tion, the persistence length of which is larger for fragile 
than for strong glass formers For non-strong liquids, 
the existence of spatially heterogeneous dynamics (SHD) 
is well established by experiment la H , |£j and Simula- 
tion H EBJiMlJipyil El M m Further, it was 
observed nil UJ, UA Mjl that groups of particles follow 
one another along string-like paths. GC proposed Q 
that dynamic facilitation and persistence of particle-flow 
direction manifest themselves in this string-like motion. 
For strong liquids, like silica, experimental studies are 
rare due to the high T g of these network-forming mate- 
rials. However, it is usually argued that the structural 
relaxation can be understood as simple activated bond 
breaking process. 

Here, we perform molecular dynamics (MD) simula- 
tions to study the relaxation dynamics of viscous silica, 
the prototype of a strong glass former In particu- 
lar, we critically test the ideas of the GC model. For 
our simulations, we apply the BKS potential 18], which 



is commonly used to reproduce structural and dynam- 
ical properties of viscous silica fl9l Eol Ell El El Eif. 
We show that the relaxation dynamics in BKS silica are 
spatially heterogeneous, indicating a complex nature of 
this dynamical process. Further, we demonstrate quanti- 
tatively that dynamic facilitation is not only relevant in 
spin models [25j, but also in a viscous equilibrium liquid. 
Finally, we show that string-like motion is less important 
in BKS silica than in non-strong liquids, suggesting a re- 
duced persistence of particle-flow direction, in agreement 
with the GC model. To arrive at these results, we follow 
Horbach and Kob [2IE3] ] and perform simulations in the 
NVE ensemble with p = 2.37g/cm 3 and TV = 8016, which 
is sufficiently large so as to avoid finite size effects at the 
studied temperatures. Further simulation details can be 
found in Refs. El El El . 



Consistent with experimental results |27j, the temper- 
ature dependence of the transport coefficients in BKS 
silica changes from a non-Arrhenius behavior at high T 
to an Arrhenius behavior at low T, where the crossover 
occurs in the vicinity of Tmct This result 

is confirmed by our analysis, see Fig. UTa). Using the 
literature value Tmct = 3330 K 22, 23 , we find that the 
diffusion coefficients D for both oxygen and silicon vary 
according to a power-law 1], D(x(T-T MC t)\ at high 
T, whereas Arrhenius laws are observed at T < Tmct- 
The activation energies E a = 4.7eV and E a = 5.0eV for 
the oxygen and silicon atoms, respectively, are consistent 
with previous results from simulations fil El | and ex- 
periments [2^, |2|| . Following the time evolution of the 
mean square displacement, r 2 (At), in Fig.^b), three dis- 
tinct time regimes can be distinguished at sufficiently low 
T [22|. Ballistic and diffusive motion at short and long 
times, respectively, are separated by a pronounced caging 
regime [lj, where the particles are temporarily trapped 
in cages formed by their neighbors. 

To demonstrate that dynamics in silica are spatially 
heterogenous, we show that highly mobile particles form 
clusters larger than predicted from random statistics. 
Following previous work [ll HI 111 flfij ] , we characterize 
the particle mobility in a time interval At by the scalar 
displacement and select the 5% most mobile particles for 
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FIG. 1: (a) Diffusion coefficients D(l/T) for the oxygen and 
silicon atoms. Dashed lines: Power- laws D oc (T — 3330 K) 7 
with 7 = 2.0 and 7 = 2.1 for the oxygen and silicon atoms, 
respectively. Solid lines: Arrhenius laws with activation en- 
ergies E a = 4.7 eV for oxygen and E a = 5.0 eV for silicon, 
(b) Mean square displacement, r 2 (At), for the silicon atoms 
at T = 3230 K. (c) Normalized mean cluster size Sw(At), 
mean string length Lw(At) and the dynamic facilitation ef- 
fect F(At), cf. Eq.H for the silicon atoms at T = 3230K. 



further analysis. Then, we define a cluster as a group of 
the most mobile particles that reside in the first neigh- 
bor shell of each other. From the probability distribution 
Ps{n; At) of finding a cluster of size n within a time in- 
terval At we calculate the weight-averaged mean cluster 
size 



Sw (At) 



E„ n*P s {n-At) 
J2 n nP s (n;At) ' 



(1) 



This quantity measures the average size of a cluster to 
which a mobile particle belongs. To eliminate the random 
contribution, we discuss the normalized size S'w(At) = 
Sw(At)/Syy, where S^y is the value that results when 
5% of the particles are selected irrespective of their mo- 
bility. Our conclusions are unchanged, when the fraction 
of mobile particles is varied over a meaningful range |2( : il| . 

Figure Hfc) shows SV(At) for the silicon atoms at 
T = 3230 K. We see that the mean cluster size exhibits 
two maxima, indicating the existence of two spatially het- 
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FIG. 2: Normalized mean cluster size, Sw(At), for the oxy- 
gen atoms at various temperatures (5250 K, 4330 K, 3870 K, 
3690 K, 3390 K, 3230 K, 3030 K). 



erogeneous dynamical processes. The secondary and pri- 
mary maxima are located at the crossovers from the bal- 
listic to the caging regime and from the caging to the dif- 
fusive regime, respectively. In Fig.0 we depict the tem- 
perature dependence of S'vy(At) for the oxygen atoms. 
The behavior of the primary maximum is an alog ous to 
that for simple liquids and polymer melts 0, llfij . Both 
the peak time ts and the peak height S\y(ts) increase 
upon cooling. Here, the growth of the clusters continues 
in the temperature range where D follows an Arrhenius 
law, cf. Fig.^a). For T«Tmct, the mean cluster sizes 
Sw(ts) for the oxygen and silicon are located at the lower 
end of the spectrum of values reported for non-strong 

liquids in M M, m 

Hence, at times when the par- 
ticles escape from their cages, BKS silica exhibits SHD, 
which is at least qualitatively comparable to that of non- 
strong liquids. The position of the secondary maximum 
is temperature independent so that both peaks merge 
for T> 4000 K. In studies of water 16], an analogous 
secondary peak was observed and ascribed to "strong 
correlations in the vibrational motion of first-neighbor 
molecules" . For BKS silica, a Boson peak was found at 
the crossover from the ballistic to the caging regime [2l| • 
Therefore, we speculate that this phenomenon leads to 
the secondary peak and, hence, can be characterized as 
localized, cooperative vibrational motion. 

String-like motion was found to be an important chan- 
nel of relaxation in fragile liquids 

mill mm t ° 

study this type of motion for silica, we follow Donati et 
al. 1 1 1 1 ] and construct strings by connecting any two par- 
ticles i and j of the same atomic species if 

min [\fi(t )- fj (t +At)\,\n (t + At) - fj {t )\]<5. 

is valid for the particle positions at two different times. 
To ensure that one particle has moved and another par- 
ticle has occupied its position, S is set to about 55% of 
the respective interparticle distance. In Fig.^c), we in- 
clude the weight-averaged mean string length Lw(At), 
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FIG. 3: Probability distributions Pm(v, At) of finding an oxy- 
gen atom v times mobile when dividing r a into time intervals 
At and identifying the mobile particles in each interval. The 
dashed lines are the corresponding expectations, Pm(v, At), 
for a random selection of mobile particles in each interval. At 
the studied temperature T' = 3030K, we find is ~10ps. 



calculated for the silicon atoms at T = 3230 K in anal- 
ogy with Eq. ^ From the small values Lw(At) « 1, it 
is clear that string-like motion is of little relevance for 
the relaxation of the silicon atoms. Though the strings 
are somewhat longer for the oxygen atoms, e.g. a max- 
imum value Lw(tL) = 1.5 is found at T = 3230 K, they 
are still much smaller than in non-strong liquids, where 
L w (t L ) > 3 for T w T M ct HI E3- Thus > we conclude 
that string-like motion in BKS silica is suppressed due to 
the presence of covalent bonds at the studied T. 

Having established the existence of dynamic hetero- 
geneities, we now measure their lifetime. For this pur- 
pose, we identify the most mobile particles in back-to- 
back time intervals, At 12 = At and At23 = At, and cal- 
culate the probability P(At) of finding a mobile particle 
in Ati2 still mobile in At23. For both atomic species at 
T = 3030K, we find a nearly constant value .P(At)«0.15 
for At in the caging regime, followed by a decrease to- 
wards the long-time limit P(At— >oo) = 0.05. These small 
values indicate that the fluctuations of the dynamical 
state are rapid. For the dynamic heterogeneities in non- 
strong liquids, such short lifetimes were observed both in 
simulation studies near Tmct fl3 1 and in experimental 
work near, but not too close to T g 

Next, we investigate the number of times a particular 
particle is mobile on the time scale of the a-relaxation 
time t q . For this analysis, we divide r Q into several 
time intervals At, identify the most mobile particles in 
each interval, and determine the probability distribution 
Pm(v, At) of finding a particle mobile in v not necessarily 
consecutive intervals. Figure[3]shows this distribution for 
the oxygen atoms at T = 3030 K. For comparison, we in- 
clude the distribution P^j(u, At), resulting for a random 
selection of 5% of the particles in the intervals At. For 
At>t$, we find Pm (v) ~ P^ { v ) hence, the subsets 
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FIG. 4: Probability distribution Pf{t, At) of finding a small- 
est distance r between an oxygen atom, which becomes mo- 
bile in At23 — At, and any oxygen atom, which was mobile in 
Afi2 = At (T = 3030 K). Pp(r,At) results when the mobile 
particles in At23 are randomly selected from the non-mobile 
particles in Ati2. The distance is calculated based on the con- 
figuration at the time t2, separating both time intervals. The 
inset shows the temperature dependence of -F(At) quantifying 
the strength of dynamic facilitation for the oxygen atoms, cf. 
Eq.H The temperatures are 5250 K, 4330 K, 3870 K, 3690 K, 
3390 K, 3230 K and 3030 K. 



of mobile particles in different intervals are statistically 
independent. Thus, high particle mobility is not limited 
to certain regions of the sample, but all particles can be- 
long to a cluster of mobile particles on the time scale 
of r a . Some deviations between Pm{v) and P^iv) for 
At<ts suggest that the relaxation of a particular region 
of the sample is not complete on this short timescale, 
consistent with the continuing growth of the clusters. 

Finally, we study whether mobility propagates contin- 
uously, as proposed in the GC model, or whether it de- 
velops at random positions, as suggested by Stillingcr 
and Hodgdon [3l| in the model of flickering fluidized do- 
mains. To answer this question, we analyze the relative 
positions of particles that are mobile in back-to-back time 
intervals At\2 = At and At23 = At, respectively. Specifi- 
cally, we calculate the probability distribution Pp (r, At) 
of finding a smallest distance r between a particle that is 
mobile in At23, but not in Ati2, and any of the mobile 
particles in Ati2. In Fig. this distribution is shown 
for the oxygen atoms at T = 3030K, where Atwtg. For 
further analysis, we compare these data with the distri- 
bution Pp(r,At), which results when the particles used 
for analysis are randomly selected from the non-mobile 
oxygen atoms in Ati2 . We see that the probability that a 
neighbor of a mobile oxygen atom becomes mobile in the 
subsequent time interval is significantly enhanced. Simi- 
lar results are observed for the silicon atoms and, hence, 
we can conclude that mobile particles in BKS silica assist 
their neighbors to become mobile, demonstrating that 
dynamic facilitation occurs. 
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The strength of this effect can be quantified based on 
the quantity 



F(At) 



f^P F (r,At)dr 
f^P*(r,At)dr' 



(2) 



where r m i n is the first minimum of the pair correlation 
function. This quantity measures by how much the prob- 
ability that a neighbor of a mobile particle becomes mo- 
bile is enhanced, as compared to the value for a random 
development of particle mobility. Figure ^c) displays 
F(At) for the silicon atoms at T = 3230 K. We see that 
the effect of dynamic facilitation is strongest at a time 
tp ss ts, i.e., when the clusters of mobile particles are 
largest. Further, F(t^)«2 indicates that dynamic facil- 
itation plays a significant role. The inset of Figure 0] 
depicts the temperature dependence of F(At) for the 
oxygen atoms. It is clearly seen that the relevance of dy- 
namic facilitation strongly increases when T is decreased, 
suggesting that this effect becomes very important near 

We have shown that the most mobile particles in BKS 
silica form clusters that are biggest at the crossover 
from the caging regime to the diffusive regime, indicat- 
ing that this strong glass former exhibits SHD. Since 
the most mobile particles show displacements compara- 
ble to the interparticle distance at these times |2f|, we 
conclude that the structural relaxation of the studied 
model liquid is not a result of a simple bond breaking 
process occurring at random positions. Further, we have 
found that the time and temperature dependence of the 
mean cluster size are similar to the behavior for sim- 
ple liquids 0, Q, 0] , implying that strong and fragile 
glass formes exhibit no principal, but at most quanti- 
tative differences in the nature of SHD on intermediate 
timescales. However, cooperative string-like motion is 
less relevant in BKS silica than in non-network form- 
ing liquids 0, 0, 0], suggesting that this dynami- 
cal pattern is suppressed by the presence of covalent 
bonds. Additionally, we have shown that the dynamic 
heterogeneities in BKS silica are short-lived and spread 
throughout the sample on the time scale of r Q . In doing 
so, the probability that a particle becomes mobile is en- 
hanced in the vicinity of another mobile particle. This 
tendency for continuous propagation of mobility demon- 
strates the relevance of dynamic facilitation. The effect 
of dynamic facilitation is strongest at times when the 
mean cluster size is maximum and it increases upon cool- 
ing, suggesting that this effect is of particular importance 
when T g is approached. 

In conclusion, we have shown for a strong glass former 
that (i) dynamics at intermediate timescales are spatially 
heterogeneous, (ii) dynamic facilitation is important and 
(iii) the persistence length of particle-flow direction as 
measured by the relevance of string-like motion is smaller 
than in fragile liquids. These points are central ideas of 



the GC model 0,0 and, hence, our findings support this 
new approach to the glass transition phenomenon. In 
addition, our results provide insights into the relaxation 
dynamics of one of the most important glass formers, 
which arc difficult to obtain from experimental studies. 

M. V. acknowledges funding of the DFG through the 
Emmy Noether-Programm. 
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